%% This file computes some of the statistics requested from the transport
% group at the world bank in April 2021

% for each country, check the share of differentiated producers that
% are along the coast and could potentially be ports 
folders;

% make a map of the consumption gains across all countries
country_list_short;

gamma = 0.10;
beta = 0.13;

row = 1;
for country_n= 1:Ncountries
    
    [ country_n ]
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name
    
    x_ver_hor = get_threshold( country );
    x_diag  = x_ver_hor;
       
    % first load in the country's grid
    load([path_load_grids, country, '_grid_', num2str(x_diag), '_', num2str(x_ver_hor), '_0.5_connected1_LAC_BASE.mat'])        
    discretized_roads=country_graph.discretized_roads; 
    unique_edges=country_graph.unique_edges;
    
    % second load in the counterfactual results for expansion 
    load([path_save_cfactuals, country, '_diag', num2str(x_diag), '_hor', num2str(x_ver_hor), '_a1_rho0_alpha0.4_sigma5_beta0.13_gamma0.1_nu1_mobil0_cong1_ngoods10_cellsize0.5_osm0_WP0_cfactual_exp_eng_BASE.mat'])
    g = cfactual.g;
        
    % from counterfactual compute change in consumption by grid cell under
    % expansion
    N = length(country_graph.places_grid);
    x = ones(N,1);
    
    for i=1:N 
        country_graph.gridmap(i).chat = cfactual.results_cfactual.Cj(i)/cfactual.results_actual.Cj(i);
    end
    
    % lanes in calibrated and counterfactual model
    I_obs = zeros( length(unique_edges),1 );
    I_exp = zeros( length(unique_edges),1 );
    for j=1:length(unique_edges)
        [j]
        I_obs(j) = g.avI( unique_edges(j,1),unique_edges(j,2) );
        I_exp(j) = cfactual.results_cfactual.Ijk( unique_edges(j,1),unique_edges(j,2) );
    end
                        
    % get discretized roads 
    for k=1:length(discretized_roads)
     %   discretized_roads(k).dI = I_exp(k) - I_obs(k);
        discretized_roads(k).Ihat = I_exp(k)/I_obs(k);
    end
    
    % Now construct tau hat based on the ratio of Ihat and Qhat
    % so we'll have 11 observations, one for each good + 1
    N = 11;
    tauhat = ones(length(unique_edges), N);    
    for n=1:N
        Q_obs = ones(length(unique_edges), 1);
        Q_exp = ones(length(unique_edges), 1);
        for j=1:length(unique_edges)
            Q_obs(j) = cfactual.results_actual.Qjkn(unique_edges(j,1),unique_edges(j,2), n);
            Q_exp(j) = cfactual.results_cfactual.Qjkn(unique_edges(j,1),unique_edges(j,2), n);
            I_obs(j) = g.avI( unique_edges(j,1),unique_edges(j,2) );
            I_exp(j) = cfactual.results_cfactual.Ijk( unique_edges(j,1),unique_edges(j,2) );
            tauhat(j, n) = ((I_exp(j)/I_obs(j))^gamma)/((Q_exp(j)/Q_obs(j))^beta);
        end
    end  
    
    discretized_roads = rmfield(discretized_roads,'opt_path');
    % shapewrite(discretized_roads, [path_raw_mapdata, '/gridmaps/', country, '_roads.shp'])  
    %shapewrite(discretized_roads, [path_raw_mapdata, '/gridmaps/', country, '_roads.shp'])  
    
    clear s
    for k=row:length(country_graph.gridmap)
        s(k).X = country_graph.gridmap(k).X;
        s(k).Y = country_graph.gridmap(k).Y;
        s(k).chat = country_graph.gridmap(k).chat;
        s(k).Geometry=country_graph.gridmap(k).Geometry;
    end
    
    % increment the counter for the next country 
    
    % shapewrite(s, [path_raw_mapdata, '/gridmaps/', country, '.shp'])  
end 

shapewrite(discretized_roads, [path_raw_mapdata, '/gridmaps/', country, '_roads.shp']) 
